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Abstract 

The biological effects of interventions to control infectious diseases typically depend on the intensity of pathogen 
challenge. As much as the levels of natural pathogen circulation vary over time and geographical location, the development 
of invariant efficacy measures is of major importance, even if only indirectly inferrable. Here a method is introduced to 
assess host susceptibility to pathogens, and applied to a detailed dataset generated by challenging groups of insect hosts 
(Drosophila melanogaster) with a range of pathogen {Drosophila C Virus) doses and recording survival over time. The 
experiment was replicated for flies carrying the Wolbachia symbiont, which is known to reduce host susceptibility to viral 
infections. The entire dataset is fitted by a novel quantitative framework that significantly extends classical methods for 
microbial risk assessment and provides accurate distributions of symbiont-induced protection. More generally, our data- 
driven modeling procedure provides novel insights for study design and analyses to assess interventions. 
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Introduction 

Hosts exposed to disease-causing agents respond in accordance 
to the challenge dose. Therefore dose-response curves contain 
information about disease processes that can be extracted by 
suitable analytic frameworks. Early examples concerning micro- 
bial risk assessment include counting lesions caused by tobacco 
mosaic virus on plant leaves [1], as well as human responders to 
experimental challenge with polio viruses [2], Vibrio cholerae [3] 
and Streptococcus pneumoniae [4], for escalating challenge doses. 
Dose-response models have been in use for analyses and 
extrapolation of experimental datasets [5] . 

Models that account for the sigmoidal shape in log-linear scale 
of the typical dose-response curve have been derived mechanis- 
tically, based on the assumption that each individual pathogen has 
a probability of infection independent of others, the so-called 
independent action hypothesis [6] . This results in a one-parameter 
exponential-function model [7]. The frequent observation of 
shallower-than-exponential, or overdispersed, relationships has 
then prompted the implementation of heterogeneity in the 
probability of infection of individual hosts [8-10]. 

In the 1960s, Fummoto and Mickey [9] developed a dose-response 
model that could accommodate both shallow and steep increases in 
the response by considering the probability of infection of individual 
hosts described by a Beta-distribution. Although a mechanistic 
justification for this specific distribution has not been given, the model 
has been widely applied in microbial risk assessment due to its ability 
to outperform the simple exponential model [5] . 

Susceptibility distributions other than Beta have also been 
considered and are more commonly used in frailty models adopted 



in survival analysis [1 1], where the data consist of survivor counts 
over time in host groups that are constantly subject to a hazard 
[12,13]. These frailty models appeared in the 1980s and have since 
been adapted to infection hazards, where surviving signifies 
remaining uninfected [14-16]. While most informative when the 
exposure is continued or repeated over time, these formalisms 
would be inadequate for estimating distributions of susceptibility to 
infection from instantaneous challenge protocols. 

The importance of accounting for time between challenge and 
observable toxicity responses to pathogens or other agents has 
been recognized. Recent models in ecotoxicology [17,18], 
consider explicit kinetics within exposed organisms. Also in 
microbial risk analysis, previous studies [19,20] have included 
time postinoculation as an additional parameter in classic dose- 
response models, although using an approach that conceptually 
allows for a different susceptibility distribution at each time point. 
Here we present a schema to infer a distribution of host 
susceptibilities to infection that holds consistently across dose 
and time. We introduce an experimental design and inference 
framework that enables such inferences by analyzing simulta- 
neously a collection of survival curves, each representing a 
different challenge dose. The resulting Beta distributions are 
compared against those obtained by classic dose-response models 
based on single day measurements. 

Recent evidence for symbiotic interactions that reduce host 
susceptibility to pathogens has stimulated the development of 
quantitative frameworks to assess the levels of individual and 
population protection attributable to specific symbionts. The 
intracellular bacterium Wolbachia, found among many arthropod 
species including Drosophila melanogaster, is one such symbiont 
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Author Summary 

While control options for plant, animal, and human 
pathogens are emerging rapidly, reliable assessment of 
the effect of interventions in biological systems presents 
many challenges. A major question is how to connect 
laboratory experiments and measurements with the 
relevant process in natural settings, where hosts are 
subject to pathogen exposures that vary in time and 
geographical location. With this aim, measures of protec- 
tion that are invariant under varying exposure intensity 
need to be developed and integrated with mathematical 
models. In this article, we introduce a method to assess 
host susceptibility to pathogens, and apply it to survival of 
Drosophila melanogaster challenged with different doses 
of Drosophila C virus. By replicating the procedure in 
groups of flies that carry the symbiont Wolbachia, we are 
able to estimate how the viral protection induced by this 
intracellular bacterium is distributed in the host popula- 
tion. Our results disentangle host infection status from 
observed mortality, accounting naturally for time since 
exposure. The multiple-dose design proposed challenges 
traditional study designs to assess interventions. 



[21,22]. To analyze the protection conferred by Wolbachia to D. 
melanogaster, we apply our inference framework simultaneously to 
two sets of time-dependent dose-response data: in one set the flies 
carry the symbiont bacterium Wolbachia (Wolb + ); while in the 
other they do not (Wolb~). In this instance we extract the Beta 
distribution that best describes individual protection attributable to 
Wolbachia, as well as population statistics valid across entire dose 
ranges. 

Methods 

Survival data 

We used virus free D. melanogaster lines with DrosDel w 1118 
background, with or without the endogenous Wolbachia strain 
wMelCS [21,23,24]. Flies were reared in standard food at 25°C. 
To assure that potential for heterogeneities are minimized by the 
experimental procedure, we used fifty 3-6 days old adult males per 
group, 10 per replicate and 5 replicates. To study the response to 
viral infection, we anesthetized with C0 2 and pricked flies with 
different doses of Drosophila C virus (DCV). We used tenfold serial 
dilutions - from 10 10 TCID 50 /ml to 10 4 TCID 50 /ml - in Tris- 
HC1 buffer, pH 7.5. Controls were pricked with buffer solution 
only. We used the pricking protocol described in [24], produced 
and titrated virus as in [21]. After pricking, we kept flies at 18°C 
and checked daily survival until day 80 and twice a week until the 
end of the experiment. Food was changed every 5 days. We 
summarized the data in 16 dose-response curves (8 per group, 
including control) from day 0 after treatment until day 139 
(Dataset SI). 

Dose-response model 

Starting from established models, we refine the occurrence of 
mortality from infection, i.e. the response, as a function of the 
concentration of infectious units given to hosts, i.e. the dose. We 
present a step-by-step derivation of descriptions that integrate 
dimensions that are usually treated separately as well as the 
motivations for doing so. 

Assuming independent action of infectious units, each unit 
has probability p of causing an infection, while for d infectious 
units infection occurs with a probability described by 



Binomial(d, p). Given further considerations about the distri- 
bution of infectious units in a homogeneous solution (see [9] for 
a complete derivation of the expression), the number of units 
causing infection can be described by a Poisson distribution, 
resulting in the exponential dose-response model [7], that 
describes the probability of infection in a host challenged with 
pathogen dose d: 



Whom = 1 — e 



-pd 



(1) 



This most basic formulation is hereafter referred to as the 
homogeneous dose-response model. 

Furumoto and Mickey [9] expanded this formulation by 
allowing the probability of infection to be described by a 
parametric distribution, specifically the Beta distribution. To 
facilitate normalization across datasets, here we maintain the 
probability /; fixed across individual hosts (as in [25]), and 
introduce a multiplicative parameter, the susceptibility factor 
0 < x < 1 , to describe any natural or induced effect that decreases 
susceptibility. We assume that susceptibility to infection is Beta- 
distributed so as to describe the variation of susceptibility in the 
host population. Thus, we obtain the probability 7thet that a host 
contracts infection as 



Whet = 1 ' 



xpd q(x)dx, 



(2) 



where q{x) - 



' B{a,b) and B is tire Beta function. We 



refer to this formulation as the heterogeneous dose-response model. 

At last we introduce a small parameter £ to account for a small 
probability of ineffective challenge, such that M ~ Binomial (n, 
( 1 — e) n) is the random variable representing the number of 
infected hosts, in a group of n hosts challenged with a given dose. 
Assuming that an ineffectively challenged host behaves like a control 
host with regard to death rates, the probability that m hosts are dead 
a number of days after challenge is then 



P(M=m)-- 



[(l-e)n] m [l-(l-e)n}' 



(3) 



where n is either 7ihom (1) or 7ih e t (2) depending on which dose- 
response model is adopted. 

The parameters to be estimated for this dose-response model 
are the maximum probability of infection per infectious unit (p), 
the shape parameters for the Beta distribution that describes the 
susceptibility factor (a, b), and the probability of ineffective 
challenge (e ). 

These models require a choice of how many days post-challenge 
cumulative mortality should be measured, which is difficult to 
establish for host-pathogen systems where times to death from 
infection or other causes overlap significantly. To overcome this 
difficulty, we develop a model that integrates an explicit 
representation of time to death with the dose-response process 
for infection just described. It should, however, be noted that time 
is introduced with the main purpose of enabling the use of survival 
curves to obtain robust estimates for probabilities of infection given 
different challenge intensities and consistently infer susceptibility to 
infection. From this perspective, parameters defined from now on 
should be regarded as auxiliary and will be implemented as simply 
as possible. 
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Time-dependent model for control group 

We first consider a survival model for a control group of flies 
pricked with buffer solution only (no Z)CV), subject to two 
hazards: ho, an age-dependent death hazard rate; and h K , a 
background age-independent death hazard rate. The overall death 
hazard rate for uninfected hosts is therefore 



h v (t) = h 0 (t) + h, c (t). 



(4) 



Denoting Tjj the random variable representing time to death of 
control hosts, we have 



than deaths due to aging, we constrain the mean time to infection 
death to be lower than old-age death, i.e. m\<mo, and the 
probability of dying before the end of the study to be greater for 
infected hosts, i.e. P(Tq < t m . dx ) <P(T\ < ? milx ), where ? max is the 
last day of the experiment. 

To construct the likelihood to be maximized by the parameter 
estimation procedure, we let Dj be the random variable denoting 
the day fly j died and S the random number of survivors up to 
'max- Then the likelihood of observing the actual number of 
survivors 5 and the times of death d=[d\, ... , d„ _ s ], for a 
given dose is 



P(T u = t)=P(T K = t)P(T 0 >t)+P(T 0 = t)P(T K >t), (5) P( S = s,D = d) = Y j P(S = s,D = d\I = i)P(I = i) 



where To and T K are the times to death from ho and h K , 
respectively. Their corresponding distributions are assumed to be 
To~Gamma(mo, so) and T K ~ Uniform(\/K) , where k is the 
background mortality rate, mo is the mean time to death, and So is 
the shape parameter for the Gamma distribution of day of death 
from aging. 

Time-dependent dose-response model 

Hosts challenged with pathogen can become infected or remain 
uninfected and this infection status is hidden. If uninfected, they 
are subject to the age-dependent hazard rate that affects control 
hosts, h\j; if infected, they are subject to an infection hazard rate, 
hi, and the age-independent background mortality. Thus the 
overall hazard rate of infected hosts is 



h I (t) = h l (t) + h K (t). 



(6) 



Now let / ~ Binomial(n, (1— e)ti) be the random variable 
representing the number of hosts infected by challenge with a 
given pathogen dose. Then the probability that i hosts are infected 
after n hosts were challenged is 



/»(/ = /)= ^ J [(l- 8 )7 I ] , [l-(l-8)7 [ r«, (7) 

where n is either 7ihom (1) or 7ihet (2) depending on which dose- 
response model is adopted. 

Let T be the random variable representing the time to death of 
hosts challenged by a given pathogen dose. The probability density 
of observing a death event at time t given that i hosts are infected is 

P ( T =t\i=i) = i ^ l) -P(T u = t) + -P(T I = t), (8) 

n n 

where Tj denotes the distribution of time to death of infected 
hosts, given by 

P{ T, = t) = P( T K = t) P( r, > t) + P( Ti = t) P{ T K >t) (9) 

and T\ is the distribution of times to death from the infection 
hazard rate h\ . This distribution is assumed to follow 
T\ ~ Gamma{m\, s{), where m\ is the mean time to death of 
infected hosts, and S\ is the shape parameter for the Gamma 
distribution of day of death from infection. 

In setting the priors for parameter estimation we note that 
background mortality is small and therefore K is kept small by 
setting 1/k to be much greater than the last day of the 
experiment. To enforce that deaths due to infection occur earlier 



E 



P(T > fmax|/ = if 11 P(dj -\<T<dj\I = i) 
7 = 1 



(10) 

P(l=i) 



Since the observations for each dose are independent, taking the 
product of the likelihoods over the different doses yields the global 
expression for the likelihood of the entire dataset. 

In this time-dependent dose-response model, the parameters to 
be estimated are the maximum probability of infection per 
infectious unit (p) used for normalization purposes, the Beta 
distribution shape parameters to describe variation in susceptibility 
factor (a,b), the parameters that control death due to aging (mo, so), 
infection (mi, Si), and background mortality (k ), as well as 
probability of ineffective challenge (e ). Parameters K and e are 
typically small and were introduced to improve performance of the 
likelihood. 

Parameter estimation 

Model parameters were estimated using Markov chain Monte 
Carlo sampling implemented with the PyMC package [26] (code 
available from [27]). The prior distributions considered are 
listed in Table 1 . Initial values were chosen so as to start with a 
non-zero likelihood. Using Metropolis-Hastings algorithm, we 
ran two separate chains for 252,000 iterations. The first 27,000 
iterations were discarded. The recording interval was set to 250 
so that the autocorrelation between samples was negligible. 
Convergence was assessed by inspection of the trace plots. All 
analyses were performed on the pooled samples from the two 
replicate chains. 

Results 

Groups of Wolbachia-negative (Wolb~) and positive (Wolb*) D. 
melanogaster flies were challenged with a range of DCV doses and 
survival curves were traced as shown in Figure 1 . This dataset was 
analyzed by applying the models introduced in Methods. 

Susceptibility distribution from selected day mortality 

To emphasize the importance of day selection to infer 
distributions of susceptibility to infection by classic dose-response 
models [5] we have applied these procedures to mortality data 
observed by two specific days (30 and 50). Parameter estimates 
from these models are listed in Table 2. The model fits to the 
mortality data at the selected days are shown in Figure 2, as well as 
the associated distribution of Wolb + susceptibilities and the 
posterior samples for the Beta distribution shape parameters. For 
simplicity we have adopted the homogeneous model for Wolb 
and focus on comparing susceptibility distributions of Wolb + 
inferred at different days. Mean protection conferred by 
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Table 1. Model parameters and their corresponding prior distributions. 



Symbol 


Meaning 


Prior 


m 0 


Mean time to death from aging 


U(0, 140) 


Sq 


Shape of the Gamma distribution for death from aging 


U(0, 500) 


mf , mf 


Mean time to death from infection (for Wolb~ and Wolb + , respectively) 


1/(0, m 0 ) 


if, if 


Shape of the Gamma distribution for death from infection (for Wolb~ and Wolb + ) 


(7(0, 100) 


P 


Per viral particle probability of causing infection 


1/(0,1) 


a, b 


Shape parameters of the Beta distribution for the susceptibility to infection of Wolb + 


(7(0.1,10) 


K 


Background mortality rate, from causes other infection or aging 


u(io- 6 , io- 2 ) 


8 


Probability of ineffective challenge 


JV(0.001, 0.00125)[0,1] 



U[x,y) is a Uniform distribution from x to y. N{x,y)[w,z] is a normal distribution with mean x and standard deviation y truncated so its values are always between w and z. 
doi:1 0.1 371 /joumal.pcbi.l 003773.T.001 



Wolbachia in this illustration is estimated as 79% and 56%, based 
on mortality measurements at day 30 and 50, respectively. 
Moreover, the distributions have fundamentally different shapes, 
with the appearance of a high susceptibility group as time 
progresses. This sensitivity to the day by which mortality data are 
collected is a concern that raises the need to disentangle infection 
status from the associated time-dependent mortality. In the 
following sections, infection and mortality are estimated explicidy 
using the integrated time-dependent model described in Methods. 
The procedure is illustrated in Figure 3. 

Aging and background mortality 

Control curves from Wolb - and Wolb + flies pricked with buffer 
solution (no DCV) were compared with the Kaplan-Meier method 
using the log-rank test and no significant difference was found 
(with a p-value of 0.47). By fitting the uninfected time-dependent 
model (4—6) to the control survival curves (Figure 1) we estimated 
the parameters describing aging (mo, So) an d background (k) 
mortality (Table 3). 

Susceptibility distribution from survival curves 

For each group of flies (Wolb and Wolb*), the time-dependent 
dose-response model constructed in Methods was fitted simulta- 
neously to the entire dataset of survival curves (one for each DCV 
challenge dose), fixing across doses the disttibution of times to death 
from infection (m/, si) and aging (mo, So), while estimating the 



susceptibility parameters (p, a, b) that govern the dependence of 
response on challenge dose according to the adopted dose-response 
model. The estimated parameter values are listed in Table 4. The 
deviance information criterion (DIC) [28] favored the homogeneous 
model for the Wolb group and the heterogeneous model for Wolb + 
(Text SI). Mean time to death from infection is 9 and 14 days in the 
Wolb - and Wolb + groups, respectively. The variance in time to 
death from infection is lower for W olb , with a standard deviation of 
2 days, compared to 6 days in the Wolb + . Figure 4 compares fitted 
with observed survival curves. 

The fitted dose-response curves that result from this analysis are 
shown in Figure 5A, while the inferred distribution of Wolb + 
susceptibilities normalized by the Wolb measure is displayed in 
Figure 5B and the corresponding posterior distribution of the Beta 
shape parameters is in Figure 5C. Given the homogeneity in the 
Wolb group, the distribution of susceptibility in Wolb + provides a 
direct indication of how antiviral protection conferred by 
Wolbachia is distributed among its carriers. Typically defined as 
1 — RR, where RR is the risk reduction attributed to the 
susceptibility modifier (Wolbachia in this case), we determine the 
mean protection conferred by the symbiont to its host as 85% 
(with a 95% HPD of 60-93%). 

Comparison with selected day mortality 

To assess the best possible performance of classic methods [5] in 
the inference of susceptibility distributions (for Wolb + in the case) 



B 




0 20 40 60 80 100 120 140 
days post challenge 



20 40 60 80 100 120 140 
days post challenge 



control 

10 4 TCID 50 

10 5 TCID 50 

10 6 TCID 50 

10 7 TCID 50 

10 8 TCID 50 

10 9 TCID 



10 10 TCID 



oo 



Figure 1. Survival curves for Wolb (A) and Wolb + (B) groups of D. melanogaster. Dots represent experimental data. Dark blue curves show 

the model fit to the survival of control flies. Shaded areas represents 95% CI (credible intervals). 

doi:10.1371/journal.pcbi.1003773.g001 
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Table 2. Estimated parameters by applying dose-response models to selected day mortality. 



Mortality data 


Parameter 


Median 


95% HPD* 


30 dpc b 


P 


2.33 10~ 6 


[1.67 10~ 6 , 3.13 10~ 6 ] 




a 


0.30 


[0.21, 0.41] 




b 


1.10 


[0.29, 2.53] 




e 


1.78 10~ 3 


[4.90 10~ 4 , 3.49 10~ 3 ] 


50 dpc b 


P 


2.65 10~ 6 


[1.82 10~ 6 , 3.47 10~ 6 ] 




a 


0.34 


[0.24, 0.51] 




b 


0.42 


[0.12, 0.93] 




e 


1.83 10~ 3 


[3.60 10~ 4 , 3.32 10~ 3 ] 



a High posterior density interval. 

b Days post-challenge. 

doi:1 0.1 371 /journal.pcbi.1 003773.t002 



we must have previously reduced the set of survival curves to a set 
of effectively infected proportions - one entry per challenge dose. 
To search for a range of days in which absolute mortality might 
provide an approximate indication of infection, we compare the 
estimated proportions effectively infected by each challenge dose 
with the mortality proportion measured at each day. Using a 
normalized Euclidean distance between these two measures, a 
day-selection score is provided by the red curve in Figure 6. We 



identify day 30 as optimal and 1 7 — 46 as the interval of days in 
which the score is at least 95 % of the optimal. Reassuringly, the 
optimal day appears to coincide with the saturation of infection- 
induced mortality (see position of vertical dash-dotted gray line in 
relation to the Gamma distributions). We now recall Figure 2 and 
Table 2 for the inferences based on day 30 mortality data to 
confirm that classic dose-response models can in principle infer 
susceptibility distributions that are consistent with those obtained 







dose susceptibility, x estimated a 

Figure 2. Dose-response curves and susceptibility distributions inferred from mortality measurements 30 and 50 days post- 
challenge. Dose-responses models adopted here are the standard formulations (1-3). A,D, Curves represent the fitted dose-response model to 
mortality on selected day post-challenge (dots), for Wolb~ (black) and Wolb + (blue). Shaded areas represent the 95% CI. B,E, Distribution of 
susceptibility to infection in Wolb + . The posterior median distribution is the curve and the shaded area is the 95% CI. C,F, Posterior samples of the 
Beta-distribution shape parameters describing Wolb + susceptibility in blue. Red dot mark the median of the respective distributions. The 
homogeneous model was adopted for Wolb - . 
doi:1 0.1 371 /journal.pcbi.1 003773.g002 
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2 groups of hosts 



no intervention 



intervention 




control infection challenge 
(doseO) (doses 1...N) 




03 
> 

> 
13 



control infection challenge 
(doseO) (doses 1...N) 




time 



survival curves 




Aging and baseline mortality independent of infection 



2. Infection parameters and host susceptibility distributions 



Figure 3. Schematic illustration of the proposed experimental design and inference procedure. 

doi:10.1371/journal.pcbi.1003773.g003 



under our extended model (Figure 5). A major issue, however, is 
that results are sensitive to a day-selection criterion that relies on 
having previously carried out the entire procedure. The appearance 
of a high susceptibility group in distributions inferred at later days 
are an artifact due to the accumulation of background mortality that 
should be factored out. These results highlight the importance of 
adequately representing the time dimension in the analysis. 

Discussion 

Dose-response models have become standard quantitative 
frameworks in microbial risk assessment. Less recognized is their 
ability to estimate host trait distributions. Here we illustrate the 
concept by extracting host susceptibility distributions from 
mortality measured as a function of pathogen challenge dose, 
but similar procedures can be developed for measures of infection 



Table 3. Estimated parameters governing time to death from 


causes other than DCV infection. 






Parameter 


Median 


95% HPD 


m 0 


117.18 


[114.99, 119.84] 


■so 


118.93 


[80.19, 166.15] 


K 


1.14 1(T 3 


[5.36 1CT 4 , 1.96 1CT 3 ] 



doi:1 0.1 371 /journal.pcbi.1 003773.t003 



or infectiousness (instead of mortality), and can be made a function 
of other environmental variables such as temperature or humidity 
(instead of dose). Understanding how to detach host trait 
distributions from environmental variables is crucial for the 



Table 4. Parameters governing estimated number infected per 
dose of DCV challenge and time to death from infection using 
time-dependent dose-response models described in Methods. 




Parameter 


Median 


95% HPD 


P 


1.73 10~ 6 


[9.58 10~ 7 , 2.67 10~ 6 ] 


a 


0.47 


[0.25, 0.85] 


b 


3.21 


[0.34, 8.40] 




1.89 10~ 3 


[4.55 10~ 4 , 3.40 10~ 3 ] 


mf 


9.34 


[9.10, 9.58] 


S T 


35.79 


[26.60, 47.05] 


mf 


13.79 


[11.31, 14.94] 


>t 


5.59 


[4.70, 11.12] 


mo 


115.20 


[113.94, 116.45] 


So 


140.39 


[116.80, 166.97] 


K 


2.15 10~ 3 


[1.65 10~ 3 , 2.71 10~ 3 ] 


Parameters with superscripts and + 
respectively. 

doi:1 0.1 371 /journal.pcbi.1 003773.W04 


relate to Wolb~ and Wolb + groups, 
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days post challenge days post challenge days post challenge days post challenge 



Figure 4. Fit of time-dependent dose-response model to survival curves. Black and blue dots are the observed proportions surviving over 
time for Wolb~ and Wolb + groups, respectively. The curve is the fitted mean posterior survival over time and the shaded area is the 95% CI. Fifty flies 
per group were pricked with: A, buffer solution (shown for comparison but not used on this analysis); and B, 10 4 ; C, 10 5 ; D, 10 6 ; E, 10 7 ; F, 10 s ; G, 10 9 ; 
H, 10 10 TCID 50 DCV. 
doi:10.1371/journal.pcbi.1003773.g004 



formulation of measures that can be transported between 
laboratory and natural conditions [29,30]. 

We address this problem with an experimental design and 
inference framework that enables the estimation of distributions of 
host susceptibility to infection by analyzing simultaneously a 
collection of survival curves, each representing a different 
challenge dose (Figure 3). The procedure is illustrated on a 
specifically collected dataset where two distinct groups of hosts (D. 
melanogaster) were experimentally challenged by viruses (DCV): 
one group consists of isogenic flies where no significant variability 
in susceptibility to infection is found; and another with the same 
genetic background but now carrying the symbiont bacterium 
Wolbachia known to reduce susceptibility to DCV [21,22]. 

Our inferences indicate that Wolbachia confers on average 85% 
DCV protection to D. melanogaster under the specified laboratory 
conditions, and suggest significant variability in this effect. This 



variance in susceptibility is induced by the symbiont, since model 
selection criteria did not support heterogeneity in the susceptibility 
of flies not carrying Wolbachia. Since the Drosophila and 
Wolbachia populations used in this study are isogenic, the 
heterogeneity in susceptibility of Wolbachia-csLrr/mg flies uncov- 
ered here indicates variation in the host-microorganism interac- 
tion that lacks a genetic basis. A simple hypothesis is that variance 
in Wolbachia levels at the individual host level leads to variance in 
resistance to viruses. Although several lines of evidence support 
this hypothesis [31-34], further experiments are required to 
discriminate whether heterogeneity in resistance is directly linked 
to variance in Wolbachia levels or, alternatively, a result of another 
environmental/ physiological variance that is only expressed in the 
presence of Wolbachia. 

Previous estimates of protection were based on survival analysis 
or viral titres in a dose-specific manner [21,22,24]. To our 




Figure 5. Dose-response curves and susceptibility distributions inferred from survival curves. A, Curves represent the estimated dose- 
response relationships from fitting the model described in Methods to survival over time, for Wolb~ (black) and Wolb + (blue). Shaded areas represent 
the 95% CI. B, Distribution of susceptibility to infection in Wolb + . The posterior median distribution is the curve and the shaded area is the 95% CI. C, 
Posterior samples of the Beta-distribution shape parameters describing Wolb + susceptibility in blue. Red dot marks the median of distribution. 
doi:10.1371/journal.pcbi.1003773.g005 
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Figure 6. Selection of optimal days to collect mortality 
measurements for traditional dose-response models. The 

red line traces a score for how well mortality at any given day 
represents infection estimated by the time-dependent model (refer 
to axis on the right). The score is given by Q=l — 



ber of doses in the dataset, x,— 



(2A), where A denotes the num- 

(x ; + ) represents the proportion 
infected in the Wolb~ (Wolb + ) group subject to DCV dose /, and y~(t) 
0>y"(f)) the observed mortality proportion overtime in the Wolb~ (Wolb + ) 
group subject to DCV dose j. Gray vertical lines mark the optimal day to 
measure mortality for dose-response models (day 30, dash-dotted line) 
and the limits of the acceptable range (days 17 and 46). Dashed lines 
represent the Gamma distributions that describe old-age mortality, and 
black (blue) full curves refer to the Gamma distributions that describe 
infection-induced mortality in Wolb~ (Wolb + ) (refer to axis on the left). 
Curves are the mean posterior probabilities and shaded areas represent 
the 95% CI. 

doi:10.1371/journal.pcbi.1003773.g006 



knowledge, the experimental design and analysis presented here 
provides the first estimation of protection in way that is detached 
from challenge dose. Future developments might consider: 
estimation of alternative distributions to compare with the shapes 
suggested by the Beta family; extension of the adopted experi- 
mental design to measure responses other that mortality; and 
move towards host populations and environmental conditions that 
are closer to natural systems. 

The parameters estimated here should not be seen as isolated 
from the relevant ecological context. On the contrary, they are 
intended as a first step to inform the construction of ecological and 
epidemiological models where Wolbachia, other symbionts, or 
interventions that modify host susceptibility to infection, are 
introduced to induce desired transitions in populations. The 



introduction of Wolbachia into Aedes aegypli and other arthropod 
vectors is being considered as a promising strategy to control 
dengue and other infectious diseases of humans (see [35] and 
references therein). The inference frameworks presented can be 
readily adapted to provide accurate quantification of Wolbachia- 
induced protection and integrated in population models of public 
health importance. 

The challenge of considering the time dependence of processes 
leading to observable ecotoxicity responses has also been 
addressed in toxicology where the so-called General Unified 
Model of Survival (GUTS) has been proposed [18]. These models 
simulate the time-course of external and internal processes leading 
to toxic effects on organisms to generate an output that can be 
fitted to mortality over time. While those studies tend prioritize the 
mechanistic descriptions of the toxicokinetic and toxicodynamic 
processes that damage the organisms, we have chosen to adopt a 
phenomenological approach and focus on the inference and 
interpretation of how susceptibility to infection is distributed in a 
population. 

In epidemiological systems, the baseline transmission intensity is 
often not direcdy measurable but indirectly inferred in a model- 
based manner. Dose-response models, on the other hand, can 
account for experimentally controlled patterns of exposure 
[36,37]. Variation in host susceptibility to pathogens is one 
component of both classes of systems that mostly influences 
estimates of intervention impacts [29] . Therefore, building on the 
methods developed here furthers our potential to accurately 
evaluate the burden of infectious diseases and design effective 
interventions. 
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